import numpy as np
import matplotlib.pyplot as plt

# Load the chi^2 map, skipping the header row
# The options for the file are:
#   + 'modchi2map_nufitwoSK.csv' corresponding to the information used in Fig. 11 (left)
#   + 'modchi2map_nufitwSK.csv'  corresponding to the information used in Fig. 11 (right)
#   + 'modchi2map_icecube.csv'   corresponding to the information used in Fig. 15
data = np.loadtxt('modchi2map_nufitwoSK.csv', delimiter=',', skiprows=1)
sin2theta23, deltam31, modchi2 = data.T

chi2_critical_value = 4.61 # This is the 90% CL chi2 for 2 degrees of freedom

dmsq21 = 0.0741 # in units of 10^{-3} eV^2
deltam32 = deltam31-dmsq21 # Convert to Delta m^2_{32}

matrix = np.array(modchi2).reshape(11, 11)
plt.imshow(matrix.T,origin='lower',cmap='Blues', vmin=0, aspect='auto',
           extent=(sin2theta23[0]-.01, sin2theta23[-1]+.01, deltam32[0]-.01, deltam32[-1]+.01) )
plt.colorbar()
plt.contour(sin2theta23.reshape(11, 11).T,deltam32.reshape(11, 11).T,matrix.T,levels=[chi2_critical_value])
plt.title(r'$\chi^2$')
plt.xticks(np.linspace(sin2theta23[0],sin2theta23[-1],11))
plt.yticks(np.linspace(deltam32[0],deltam32[-1],11))
plt.xlabel(r'sin$^2(\theta_{23})$')
plt.ylabel(r'$\Delta m^2_{32}$ [10$^{-3}$ eV$^2$]')
plt.show()
